# -*- coding: utf-8 -*-
"""
Created on 03 08 下午4:55 2022 
@Author : HHQUAN
@Email : 1075960398@qq.com

"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import soundfile as sf

# 两个实际系统的频响曲线，其中一个用FIR拟合需要60阶，另一个需要10阶
omega = np.arange(0, np.pi, np.pi/100)
G1 = np.zeros(omega.shape, dtype=complex)
G2 = np.zeros(omega.shape, dtype=complex)
for i in range(omega.size):
    z = np.exp(1j*omega[i])
    G1[i] = 1/(1+0.5*z**(-1)-0.2*z**(-2)+0.2*z**(-6))
    G2[i] = 1/(1+0.5*z**(-2))

noise = np.random.randn(16000) # 噪声测试成功
wav, fs = sf.read('/home/cl/PycharmProjects/Mytest01/方位估计/man_woman16k.wav')
wav1 = wav + 0.001*np.random.randn(wav.size)
b,a = [1.0], [1.0, 0.5, -0.2, 0, 0, 0, 0.2]
b1, a1 = [1.0], [1.0, 0.0, 0.5]
noise2 = signal.lfilter(b,a, wav1, -1)
plt.figure()
plt.plot(wav1, 'b')
plt.plot(noise2, 'r')

plt.figure()
plt.specgram(noise2, NFFT=320, Fs=16000, detrend='linear', noverlap=160)

plt.figure()
plt.plot(omega/(2*np.pi), abs(G1))
plt.plot(omega/(2*np.pi), abs(G2))
plt.show()

